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Abstract 

Two-component (normal and degenerate stars) mod- 
els are the simplest realization of clusters with a mass 
spectrum because high mass stars evolve quickly into 
degenerates, while low mass stars remain on the main- 
sequence for the age of the universe. Here we ex- 
amine the evolution of isolated globular clusters us- 
ing two-component Fokker-Planck (FP) models that 
include heating by binaries formed in tidal capture 
and in three-body encounters. Three-body binary 
heating dominates and the postcollapse expansion 
is self-similar, at least in models with total mass 
M < 3 x 10 5 Mq, initial half-mass radius r/, j > 5pc, 
component mass ratio m^jm\ > 2, and number ratio 
Ni/N 2 < 300 when m 2 = 1.4 M Q . We derive scaling 
laws for p c , v c , r c , and 77, as functions of mi/m2, N, 
M, and time t from simple energy-balance arguments, 
and these agree well with the FP simulations. We 
have studied the conditions under which gravother- 
mal oscillations (GTOs) occur. If E tot and E c are 
the energies of the cluster and of the core, respec- 
tively, and t r h and t c are their relaxation times, then 
e = (E tot /t r h)/{E c /t rc ) is a good predictor of GTOs: 
all models with e > 0.01 are stable, and all but one 
with e < 0.01 oscillate. We derive a scaling law for e 
against N and mi / m 2 and compared with our numer- 
ical results. Clusters with larger rr^/mi or smaller N 
are stabler. 

Subject headings: celestial mechanics, stellar dynamics — 
globular clusters : general 
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1. INTRODUCTION 

The dynamical evolution of pre- and postcollapse 
globular clusters is influenced by many factors: the 
initial mass function, the nature and efficiency of en- 
ergy generation mechanisms, galactic tides, anisotropy 
of the velocities, the initial population of binaries, 
and stellar evolution (see, for example, Spitzer 1987). 
If the goal is to model globular clusters realistically, 
then all of these effects should be included. On the 
other hand, if the goal is a deeper theoretical under- 
standing of individual dynamical processes, then sim- 
pler models can be more instructive. 

A distribution of stellar masses affects the post- 
collapse evolution of globular clusters in ways that 
have not been fully explored. The simplest nontriv- 
ial multimass models are those with just two compo- 
nents. This simplification is drastic but not entirely 
unrealistic. Since the stellar main-sequence lifetime 
is a very steep function of stellar mass, a cluster can 
be assumed to start its dynamical evolution with a 
turnoff-point mass very similar to that currently ob- 
served. Stars with higher mass have already evolved 
into degenerate remnants (such as black holes, neu- 
tron stars or white dwarfs). Because of mass segre- 
gation (also called "mass stratification," cf. Spitzer 
1987), the inner parts of a dynamically relaxed clus- 
ter should consist primarily of the turnoff stars and 
the heaviest remnants. The black holes are probably 
much more massive than present day main-sequence 
stars, but their number is expected to be too small to 
play any important dynamical role (Kulkarni, Hut, & 
McMillan 1993). Most of the white dwarfs in a clus- 
ter are expected to have masses less than the present 
turn-off point mass. The dynamical effects of such 
white dwarf stars are expected to be unimportant. 
Kim & Lee (1997) compared the numerical results of 
two- and 11-component models, which have 3 compo- 
nents for the white dwarfs heavier than the turn-off 
point mass, and they were able to find two-component 
cluster parameters that well match the 11-component 
clusters. This implies that the presence of the white 
dwarfs heavier than the turn-off point mass does not 
harm the two-component approximation. 

Neutron stars, however, can be an important com- 
ponent in dynamical evolution of globular clusters. 
Neutron stars are about twice as massive as turnoff 
stars. The fraction of the cluster mass in the form of 
neutron stars is expected to be very small, but mass 
segregation may make neutron stars major compo- 



nent in the central parts. (For observational evidence 
supporting the dominance of heavy degenerates, see 
Gebhardt & Fischer 1995 and Phinney 1993.) As a 
first approximation, the globular clusters can be rep- 
resented by two components: neutron stars and main- 
sequence stars. 

Two-component clusters have several interesting 
features that are not present in single-component 
clusters. Mass segregation accelerates the early phases 
of core collapse significantly (e.g. Yoshizawa et al. 
1978; Inagaki & Wiyanto 1984; Spitzer 1987). The 
central parts quickly become dominated by heavy 
component, even though bulk of the mass of the clus- 
ter is in the light component. Tidal captures between 
a neutron star and a main-sequence star are frequent, 
and binaries composed of neutron stars can also be 
formed by three-body processes. These binaries will 
eventually drive the postcollapse evolution. 

The purpose of this paper is to obtain a series of so- 
lutions for the dynamical evolution of two-component 
clusters, with emphasis on the postcollapse phase. We 
study relative importance of tidal-capture and three- 
body binary formation for heating the cluster core. 
We are particularly interested in the importance of 
unequal stellar masses for gravothermal oscillations. 

Our model is quite simple compared to real clus- 
ters. Stellar evolution has not been included; it is 
most important for the very early evolution of clus- 
ters (Chernoff & Weinberg 1990, Drukier 1995). Pri- 
mordial binaries could be even more important than 
the binaries formed by dynamical processes through 
the early postcollapse phases (Gao et al. 1991), but 
we neglect them. We have ignored the tidal field of 
the galaxy, even though tidal limitation qualitatively 
alters postcollapse evolution of the cluster mass and 
radius (Henon 1961, Lee & Goodman 1995), and tidal 
shocks may further hasten the destruction of clusters 
(Kundic & Ostriker 1995; Gnedin & Ostriker 1997). 
External tidal fields are probably not important for 
gravothermal oscillations and other phenomena per- 
taining to the core, except insofar as they modify the 
total cluster mass. 

This paper is organized as follows. In we de- 
scribe the models and methods of our calculations. 
Heating mechanisms that derive the postcollapse ex- 
pansion are compared in Q In §|4|, we obtain a series 
of solutions and give simple analytic expressions for 
the evolution of parameters of postcollapse clusters. 
Gravothermal oscillation in two-component models is 
examined in da. The final section summarizes our 



results. 

2. MODELS 

2.1. Computational Method 

The dynamical evolution of collisionless stellar sys- 
tems under the influence of two-body relaxation is 
well described by the Fokker-Planck equation. We 
have restricted ourselves to isotropic models, in which 
the stellar orbital distribution function depends only 
on energy (and time). A recent study by Taka- 
hashi, Lee, & Inagaki (1997) showed that the ra- 
dial anisotropy in the halo becomes highly suppressed 
when a tidal field is imposed. The global evolu- 
tion can be simply described by an isotropic model. 
The multi-component Fokker-Planck equation can be 
written as follows: 
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where fi(E) and Fi(E) are the distribution function 
and the particle flux in energy space E, respectively, 
for the i-th component. Formally Fi can be expressed 
as 

dfi 
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where D$ and Dee are the Fokker-Planck coeffi- 
cients. The statistical weight factor p(E) is given by 



p(E) = 4 f 
Jo 



r 2 v dr, 
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where </> -1 (i?) denotes the maximum radius that a 
particle with energy E can reach in the spherical po- 
tential (j){r). 

In order to account for the effects of binaries, we 
need to modify the Fokker-Planck equation above. 
Statler, Ostriker & Cohn (1987) developed a sophis- 
ticated scheme to include the dynamical effects of 
tidal binaries for a model initially composed of a sin- 
gle component. Lee (1987) modified the method of 
Statler et al. (1987) to include both tidal binaries 
and three-body binaries. However, the number of 
dynamically-produced binaries is always a very small 
fraction of the total number of stars. The most im- 
portant effect is heating — the addition of entropy to 
the orbital distribution function — which can be sim- 
ply accounted for by modifying the particle flux in en- 
ergy space. Such an approach has been taken in many 
studies dealing with the three-body binary heating 
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(see, for example, Cohn 1985; Lee, Fahlman, & Richer 
1991). Normally the tidal binaries are more abundant 
and simple correction of the particle flux could cause 
some errors. However, we find that our approach pro- 
vide excellent agreement with more complex schemes. 

The orbit-averaged heating coefficient due to this 
heating becomes 



Hi{E) 
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where E is the heating rate per unit volume and r max 
the maximum radius accessible to a star with energy 
E. The modified particle flux then takes the form 

Bf. 

Ft (E) = - [ mi D E + Hi (E)] ft (E)-D EE -^. (5) 

The second-order coefficient (Dee) is also affected by 
binary heating, but this diffusive effect is probably 
much less important for the long-term evolution than 
the change in De- 

Heating by tidal binaries is mostly due to their 
ejection during close encounters with single stars (or 
sometimes other binaries) because the internal bind- 
ing energy of tidal binaries is much greater than the 
escape energy from the cluster. Therefore the heating 
rate per unit volume is the product of the mass in all 
three stars, the binary formation rate, and the central 
potential (Lee & Ostriker 1993): 



Etc = (m n + 2m d )a tc v re in n nd4>c 



(6) 



where n is the number density, ate the tidal-capture 
cross section, v re i the relative rms velocity between 
degenerate and normal stars, and <p c the central grav- 
itational potential. The subscript n denotes normal 
stars while d represents degenerate stars. It is as- 
sumed that the binary is ejected promptly after for- 
mation by interaction with a third star. We adopt 
the following expression for ate from Lee & Ostriker 
(1986) for tidal captures between a normal star and 
a degenerate star with md/m n = 2: 
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where w e „ is the escape velocity at the normal star's 
surface, and R n the normal star's radius. 

Binaries formed by dissipationless three-body en- 
counters ( "three-body binaries" ) are usually less tightly 



bound than tidal binaries. Therefore, in their en- 
counters with single stars, some or all of the stars 
are retained by the cluster and contribute their in- 
creased orbital energy to the distribution functions 
fi ("direct heating"). Until ejected from the cluster, 
a three-body binary releases approximately 300 kT, 
where kT is the typical kinetic energy of background 
stars. Thus the three-body binary heating rate per 
unit volume can be computed by multiplying the bi- 
nary formation rate per unit volume by 300 kT: 



E 3b = 4.21 x 10 3 G 5 



niinf 



(8) 



where the summation is over all components, and v\ is 
the mass weighted, three-dimensional, central velocity 
dispersion. The numerical coefficient has been taken 
from Cohn (1985). 

The modified Fokker-Planck equation can be ac- 
curately integrated with the numerical procedure de- 
scribed by Cohn (1979, 1980). 

2.2. Model Parameters 

Our models assume an initial cluster composed of 
main-sequence stars with mass mi and neutron stars 
with mass m.2- The total masses in the form of these 
stars are M\ and M2 respectively. The total number 
of stars is denoted by N. 

We need to specify three dimensionless parameters 
in addition to the initial density and velocity pro- 
files: the total number of stars N, the ratio TTI2 / tyi\ 
of a heavy star to a light one, and the number ra- 
tio N2/N1. The dimensional scales are determined by 
the total cluster mass M, the initial half-mass radius 
rh,i, and the stellar radius R n - The ratio m\rh,i/ MR n 
is important for tidal heating because it determines 
the ratio of the capture cross section (eq. 0) to the 
projected area of the cluster (or cluster core). The 
present study may be devided into three topics and 
each topic has its own set of models. The detailed 
model parameters of those sets are listed in Tables 
[l], ||, and |]. Note that in all our runs, the total 
mass of the heavy component, Mi = N^mi, is negli- 
gible compared to the total mass of light component, 
Mi = N\m\, and thus mi w M/Ni. The initial den- 
sity and velocity profiles are given by Plummer mod- 
els with v c i/v C 2 — 1 and p c \j p C 2 = M1/M2, where p c 
is the core density. Both three-body binary heating 
and tidal binary heating are included and clusters are 
assumed to be isolated (i.e. no tidal cut-off). 
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Our values for 7712 / ra\ range from 2 to 5. If we iden- 
tify mi with the heaviest main-sequence stars, then 
m 2 /mi = 2 is a suitable choice. Since the bulk of the 
cluster consists of stars below the turnoff, however, 
mi should be somewhat smaller than the turnoff- 
point mass. Finally, we assume a linear mass-radius 
relation such that 1 M Q corresponds to 1 i?©. 

3. HEATING MECHANISMS 

Most heating takes place in the core where the stel- 
lar and binary number densities are highest. To com- 
pare the efficiencies of heating mechanisms, we first 
derive their dependence on core parameters. From 
equations (||) and (^), E tc in the core is approxi- 
mately proportional to v n n n rid, while in a core 
dominated by degenerate stars is approximately pro- 
portional to n|v7 7 . We have assumed that 4> c oc v^ c , 
which is valid during the postcollapse phase. Even 
during the precollapse phase, the ratio of potential 
depth to central velocity dispersion varies very slowly 
compared with other quantities. Considering again 
the fact that the core is dominated by degenerate 
stars, the central heating rate by tidal binaries is pro- 
portional to the first power of the central density, and 
the central heating rate by three-body binaries is pro- 
portional to the third power. Since the central veloc- 
ity dispersion evolves much less than central density, 
the density of degenerate stars in the core is the most 
important quantity in deciding the relative efficiency 
of the two heating mechanisms. Thus the relative im- 
portance of three-body binaries increases as the core 
collapse proceeds. 

We have searched for the division of tidal and 
three-body heating dominance in the two-component 
parameter space, and its result is shown in Table [|. 
3B denotes the models whose postcollapse expansion 
is driven by three-body binaries, and TB by tidal bi- 
naries. Three-body binary heating is relatively more 
important for clusters with smaller M and N1/N2, 
and larger ru,i and m%j 'mi. Smaller N1/N2 and larger 
rail mi give larger initial n C 2/n c i which is favorable 
to three-body binary formation. On the other hand, 
a cluster with very high initial n c starts its evolu- 
tion with tidal binary heating rate high enough not 
to give a change for three-body binary heating to be 
dominant during the collapse. Kim & Lee (1997) 
found that the evolution of multi-component mod- 
els with various power-law mass functions and M 
may be realized with two-component models with 
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Fig. 1. — Heating rates by three-body binaries (solid 
line) and tidal binaries (dashed line) of run caab. The 
units of the y-axis are arbitrarily chosen. 

1.7 < to 2 < 3 and 10 < Ni/N 2 < 50 (the actual 
fraction of neutron stars in globular clusters is much 
less than 1/100, but many low-mass normal stars do 
not contribute much on the cluster's dynamical evolu- 
tion). Although the dominant heating mechanism is 
dependent on the mass function {m%lm\ and N1/N2), 
all models with typical globular cluster masses and 
sizes (M < 10 6 M o and r h ^ > 2.5 pc) and with the 
above mass function range are dominated by three- 
body binary heating in the postcollapse phase as in 
Figure |l|, which is a typical history of the heating 
rates of our models marked with 3B in Table |l|. This 
dominance of 3B over TB in most plausible parame- 
ter range confirms the conclusion of Lee (1987) that 
in multi-component clusters, three-body binary heat- 
ing is relatively more important because there are no 
tidal captures between pairs of degenerate stars. Note 
that, even for an extreme N1/N2 value of 300, three- 
body binary heating is still dominant for clusters with 
M < 3 x 10 5 Mq and rh,i > 5 pc, which are the param- 
eter ranges that the majority of the present clusters 
have. 

Although tidal binary heating is unimportant for 
postcollapse expansion of models marked with 3B, 
tidal binary heating becomes more important when 
M and N1/N2 are larger, and when rhi and rai/rai 
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Table 1 

Dominant Heating Mechanism in the Postcollapse Phase 



M 




rh,i 


(pc) 




M 




Th,i 


(PC) 




(M e ) 


1 


2.5 


5 


10 


(M e ) 


1 


2.5 


5 


10 


1 x 10 5 


m 2 


= 2, 


iVi _ 
N 2 


30 


1 x 10 5 


m 2 
m\ 


= 2, 


Ei - 

N 2 


300 


3B 


3B 


3B 


3B 


TB 


3B 


3B 


3B 


3 x 10 5 


3B 


3B 


3B 


3B 


3 x 10 5 


TB 


TB 


3B 


3B 


1 x 10 6 


TB 


3B 


3B 


3B 


1 x 10 6 


TB 


TB 


TB 


TB 


3 x 10 6 


TB 


TB 


3B 


3B 


3 x 10 6 


TB 


TB 


TB 


TB 


1 x 10 5 


>no 
m i 


= 3, 


iVl _ 
N 2 


30 


1 x 10 5 


'»2 
m l 


= 3, 


Ei - 

N 2 


300 


3B 


3B 


3B 


3B 


3B 


3B 


3B 


3B 


3 x 10 5 


3B 


3B 


3B 


3B 


3 x 10 5 


3B 


3B 


3B 


3B 


1 x 10 6 


3B 


3B 


3B 


3B 


1 x 10 6 


TB 


3B 


3B 


3B 


3 x 10 6 


3B 


3B 


3B 


3B 


3 x 10 6 


TB 


TB 


TB 


3B 



Note. — 3B is for the models whose postcollapse expansion is driven by three- 
body binary heating, and TB by tidal binary heating. 7712 = 1.4 Mq for all models. 



are smaller. For example, our two-component model 
with M = 10 6 Af Q , r h ,% = 2.5 pc, m 2 /mi = 2, 
Ni/N 2 = 300, and m 2 = 1.4 Mq is governed by 
tidal binary heating in the postcollapse phase (see 
Figure |J). It is notable, however, that the three- 
body binary heating rate increases after core collapse 
while the rate of heating by tidal binaries is decreas- 
ing. This phenomenon is not expected for single-mass 
tidal- binary-driven postcollapse clusters, in which the 
central density and velocity dispersion evolve as p c oc 
t~ XM and v c oc i" 34 (Lee & Ostriker 1993), and 
the three-body binary heating rate as plv~ 7 oc t~ - 8 . 
In a two-component postcollapse cluster, however, if 
tidal binaries dominate and equipartition holds in the 
core, then the central density decreases less rapidly 
and the central velocity more rapidly than in a single- 
component cluster (Figure ||) . In such cases, the heat- 
ing rate by three-body binaries may even increase dur- 
ing postcollapse expansion. 

The ultimate age of the models shown in Figures 
appears unrealistically long. But the evolution 
of isolated clusters slows down as expands (see 
§3.3 below), whereas r^, t r h, and M decrease during 
the postcollapse evolution of tidally limited clusters, 
which are more realistic (cf. Lee & Ostriker 1987). 
Thus, the evolutionary state of an isolated cluster at 
t, = 10 11 yr corresponds roughly to that of a tidally 



limited one at a much earlier time. Mass loss by 
overflow of the tidal boundary is much more rapid 
than ejection of stars by binaries in the core, and de- 
pletes primarily the lighter component. By decreasing 
N1/N2, this last effect further suppresses tidal bina- 
ries. 

Many of the interactions that we identify as tidal 
captures would actually have lead to mergers (Bcnz 
& Hills 1987). Ultimately, the mass of the normal 
star might still be ejected, but the neutron star would 
probably not be. Thus the parameter space for clus- 
ters whose postcollapse phase is driven by tidal bi- 
nary heating is expected to be even narrower than 
our present calculations suggest. As a limiting case, 
we have re-calculated the two-component models in 
Table |l| assuming that the all tidal capture binaries 
end up with mergers and only the masses of the nor- 
mal stars are ejected, i.e. 

Etc = m n at c v re in n nd<t>c- (9) 

The dominant heating mechanism in the postcollapse 
phase for the models in Table [j] with the above E tc 
is shown in Table |2| The dominance of tidal binary 
heating is seen only for few models, because equa- 
tion (^) gives l/(2m2/wii + 1) times less heating rate 
than equation (^) and thus three-body binaries now 
become even more important. For realistic mass func- 
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Table 2 

Dominant Heating Mechanism in the Postcollapse Phase (Merging Case) 



M 




rn,i 


(PC) 




M 






(pc) 




(M Q ) 


1 


2.5 


5 


10 


(M e ) 


1 


2.5 


5 


10 


1 x 10 5 
3 x 10 5 
1 x 10 e 
3 x 10 6 


m i 


= 2, 


ifi - 

N 2 


30 


1 x 10 5 
3 x 10 5 
1 x 10 6 
3 x 10 B 


m i 


= 2, 


JVi _ 

N 2 


300 


3B 
3B 
3B 
3B 


3B 
3B 
3B 
3B 


3B 
3B 
3B 
3B 


3B 
3B 
3B 
3B 


3B 
TB 
TB 
TB 


3B 
3B 
TB 
TB 


3B 
3B 
3B 
TB 


3B 
3B 
3B 
TB 


1 x 10 5 
3 x 10 5 
1 x 10 6 
3 x 10 6 


m i 


= 3, 


Ni _ 
N 2 — 


30 


1 x 10 5 
3 x 10 5 
1 x 10 6 
3 x 10 B 


mg 

mi 


= 3, 


JVi _ 
N 2 — 


300 


3B 
3B 
3B 
3B 


3B 
3B 
3B 
3B 


3B 
3B 
3B 
3B 


3B 
3B 
3B 
3B 


3B 
3B 
3B 
TB 


3B 
3B 
3B 
3B 


3B 
3B 
3B 
3B 


3B 
3B 
3B 
3B 



Note. — Notations and are the same as in Table |lj 



tions (N1/N2 < 50), no model in our parameter range 
(M < 3 x 10 6 and r h , ;i > pc) is marked with TB. 

4. PRE- AND POSTCOLLAPSE EVOLU- 
TION 

In this section, we show the results for 11 models 
(Group A). The parameters of these models, which 
are shown in Table ^, have been chosen to cover a 
range of plausible or instructive combinations. We 
now discuss some features of these models. 

4.1. Epoch of Corecollapse 

Isolated, single-component clusters beginning as 
Plummer models reach core collapse after a time 
t C c = 15.4 t r h,i (Cohn 1980), where t r h,i, the ini- 
tial half-mass relaxation time, does not vary much 
before core collapse. However, the ratios t cc /t r h,i 
and t cc /t rc (where t rc is the core relaxation time) 
depend strongly on the initial density and velocity 
profile (Inagaki 1985), and can be much smaller for 
choices other than the conventional Plummer model. 
Quinlan (1996) found that for single-mass clusters, t cc 
varies much less when expressed in units of t rc divided 
by a dimensionless measure of the temperature gradi- 
ent in the core. Single-mass clusters evolve by radial 
transport of energy, but energy exchange between dif- 
ferent mass components plays an important role in 



multi-component models. Mass segregation in a two- 
component model takes place initially on a timescale 
set by dynamical friction, which can be shorter than 
the corecollapse time of a single-component models 
with similar macroscopic properties. Mass segrega- 
tion and core collapse in two-component models has 
been discussed in detail by Inagaki (1985). In Table 
|3|, we have listed the t cc /t r h,i ratios for each of our 
models. 

4.2. Scaling Laws for Postcollapse Evolution 

The expansion of the core in postcollapse is de- 
termined by the dominant heating mechanism. Here 
we present scaling laws for the postcollapse evolution 
of two-component clusters driven by three-body bi- 
nary heating based on theoretical analysis and com- 
pare them with our numerical results (see Lee & Os- 
trikcr 1993 for scaling laws for evolution driven by 
tidal binary heating). 

We start with Goodman's (1993) energy balance 
analysis arguments. The energy generation rate by 
three-body binaries in the core is 



Lc w M, 



Pc 



(10) 



where the core mass M c and core radius r c are given 
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Table 3 

Parameters and Results of Simulation Group A 



Run 


mo 
mi 


Ni 

N-2 


M 
(M e ) 


A 7 


m 2 
(M e ) 




Values at t = 

Pc Vc 

(M pcr 3 ) (kins' 1 ) 


iu yr 

r c 

(PC) 


(PC) 


baab 


2 


100 


10 5 


141457 


1.4 


12.42 


1.47 x 10 5 


3.03 


0.0579 


28.8 


caab 


3 


100 


10 5 


210125 


1.4 


6.34 


1.29 x 10 5 


2.86 


0.0585 


21.2 


faab 


4 


100 


10 5 


277473 


1.4 


3.23 


1.19 x 10 5 


2.80 


0.0594 


19.9 


cdab 


3 


30 


10 5 


201299 


1.4 


3.97 


1.26 x 10 5 


2.90 


0.0600 


28.6 


cbab 


3 


300 


10 5 


212871 


1.4 


10.92 


1.26 x 10 5 


2.85 


0.0587 


21.5 


caabl 


3 


100 


10 5 


70042 


3x1.4 


6.47 


3.18 x 10 3 


2.08 


0.270 


40.9 


caab2 


3 


100 


10 5 


630374 


ixl.4 


6.28 


5.36 x 10 6 


3.91 


0.0124 


11.4 


baab3 


2 


100 


10 5 


212185 


1x1.4 


12.17 


6.10 x 10 5 


3.40 


0.0319 


22.7 


faab3 


4 


100 


10 5 


208104 


|xl.4 


3.23 


4.57 x IO 2 


2.58 


0.0884 


23.5 


caeb 


3 


100 


3 x 10 4 


63037 


1.4 


6.34 


2.21 x IO 3 


1.35 


0.210 


29.5 


cabb 


3 


100 


3 x 10 5 


630374 


1.4 


6.43 


5.61 x IO 6 


5.68 


0.0176 


15.9 



Note. — The initial half-mass radii rh,i of these runs are all 5 pc. 



by 



27r 3 

M c = —p c rl 
3w 2 



(11) 



(12) 



c 4irGp c 

On the other hand, the power required by the expan- 
sion of the cluster is 



&h ~ — y—r h . 



(13) 



Since the heavy component dominates in the core and 
the light component at rh, substitution of the core 
parameters into equation (0) yields 



L c oc mlr 3 c p 3 c v c 7 



(14) 



Slowly evolving postcollapse solutions are almost 
isothermal and the core approaches equipartition. 
Thus the Virial relation becomes 



mi 2 mi GM 



u c2 



TO 2 



m 2 r h 



(15) 



Furthermore, the temporal dependence of the above 
parameters can be obtained from the assumption that 



r h 1 
— oc . 

rh Uh 



(16) 



which would follow if the evolution were self-similar. 
Since t r h oc M 1 ^ 2 m^ 1 rff 2 , the above relation gives 



m i 



(17) 



where r/ l cc is at t = t cc , and the outer part of the 
cluster is assumed to be dominated by the light com- 
ponent. Since is almost constant until the corecol- 
lapse takes place, 



, 3/2 3/2 N Wl 

(r h ' - r - 



/' ~ 'h,i ) K ~ tcc ^ 



and for t^$>t c 



3/2 mi 



It follows from equations (16) and (|13|) that 



Eh oc m\M 



3/2^-5/2 



(18) 



(19) 



(20) 
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Fig. 2. — Heating rates by three-body binaries (solid 
line) and tidal binaries (dashed line) of our two- 
component model with M = 10 6 M Q , rh,% — 2.5pc, 
ma /mi = 2, Nx/N 2 = 300, and m 2 = 1.4 M©. The 
units of the y-axis are arbitrarily chosen. The initial 
tidal binary heating rate is so high that the three- 
body binary heating does not have a chance to take 
over. However, the three-body binary heating keeps 
increasing in the postcollapse phase. 

Demanding that the power ( pp| ) required by expansion 
be supplied by the core luminosity ([l4|), and substi- 
tuting for r h from (|l9|), one finds the following rela- 
tions for late postcollapse evolution (i.e. t 3> t cc ): 

Pc oc (n!ir 10/ Vo/3 4 -2 ; (21a) 

\mi J 

Vc oc f^r 1/2 ivi/3 M i/3 4 -i/3 ; (21b) 

\miJ 

r c oc f!!!iy /6 7V-4/3 M i/3 4 2/3. (21c) 
\miJ 

r h oc N-^Hl 1 ^^ 3 ; (21d) 

Mc oc f^V /6 7V- 2 /3M. (21c) 
\miJ 

It is notable that 77, j, one of the initial parameters of 
the cluster, is not included anywhere in equation (|2~l|). 
This, too, is a reflection of the self-similar nature of 
postcollapse evolution, which has been confirmed by 
many previous studies, but mostly for a single mass 
component. Tidally limited postcollapse evolution, 
on the other hand, is not strictly self-similar, since 



Fig. 3. — Temporal evolution of the core density 
(thick lines) and the central velocity dispersion (thin 
lines) for the same run as in Figure |[ Dotted lines 
are for the light component, dashed lines for the 
heavy component, and solid lines for the overall val- 
ues. Equipartition is approached in the postcollapse 
phase. p c is in units of M & pc~ 3 and v c in kms -1 . 

r c /rh increases with decreasing cluster mass. 

Unlike tidal binary heating, the calculation of three- 
body binary heating can be done purely dimension- 
lessly and is scalable to isolated clusters of any initial 
size, mass, and density. We can fix all of the constants 
in the scalings above from our numerical experiments: 

p c ~ 4.5xlO 5 M /pc 3 (^)" 10/3 7V 5 10 /V'^2a) 
v c ~ 3.8km/s (^.y 1 ' 2 Nl' 3 Ml'\°- 32 ; (22b) 

r c ~ 0.042 pc f^V /6 iV 5 - 4/3 M 5 1/3 *ll 8B ; (22c) 
r h ~ 35pc JV 5 _2/3 M 5 1/3 t?i 65 ; (22d) 
M c ~ 71 M Q (^y /6 N- 2/3 M 5 t^ 0B , (22c) 

where N 5 = iV/10 5 , M 5 = M/10 5 M Q , and t n = 
f/10 11 yr. Note that the exponents of t, which are 
obtained by the power-law fitting, show only small 
discrepancies from equation (]2l|). 

The postcollapse central density evolutions of some 
of our runs with the same initial conditions and the 
same N x mi /m-x values are shown in Figure ||. Ac- 
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Fig. 4. — Temporal evolution of central densities, 
three-dimensional rms central velocities, core radii 
and half-mass radii for runs baab (solid line), caab 
(dashed line), and faab (dotted line). These three 
runs have nearly the same N x m\jmi values and 
converge into the same evolutionary track in the post- 
collapse phase. p c is in units of MqPC -3 , r/j and r c 
in pc, and v c in kms -1 . 



cording to equation (21), these runs should have the 
same p c at the same £(> t cc ) although their t cc 's are 
different. After a short transition period in the be- 
ginning of the postcollapse phase, the p c 's do indeed 
converge to the same log-log slope. 

Figure ^ shows p c , v c , r c , and rh of our runs in 
Group A at t = 10 11 yr over the right-hand-sides 
in equation (plj). All panels in the figure have the 
same abscissa and ordinate scales so that a slope 
of unity represents the proportions in equation (pi]). 
As one can see from these figures, our numerical re- 
sults closely follow the scaling relationship. The best- 
fitting slopes in the figures are 1.02±0.01, 0.92±0.01, 
1.07 ± 0.01, and 0.86 ± 0.06 for p c , v c , r c , and r h , 
respectively. These slopes are close to unity as rep- 
resented by diagonal lines. The slope of rh is the 
most discrepant. This is probably due to the approx- 
imations used in equations ( flij ) through jl9|). The 
validity of equation ([l6]) will be discussed again in §||. 

According to the scalings above, the central density 
depends upon the total mass M, the epoch t, and on 



Fig. 5. — (a) Central densities, (b) three-dimensional 
rms central velocities, (c) core radii, and (d) half- 
mass radii of runs in Group A. N = N/10 5 and M5 = 
M/10 5 Mq. The straight diagonal lines represent the 
theoretical relations. 



m-2, but it is independent of Mi and M2 separately. 
Therefore, the mass function does not influence the 
central density in postcollapse (provided that there 
are enough heavy stars to fill the core — see below), 
even though it strongly influences the epoch of core 
collapse. 

The degree of concentration can be measured by 
r h/ r c which is constant during the postcollapse phase. 
From equations (27c) and (27d) we have 



r h /r c = 810 N, 



2/3 I mi 
m 2 



7/6 



(23) 



In tidally limited postcollapse models, the tidal ra- 
dius is larger than r/, by a factor of the order of 
10. This means the concentration parameter of post- 
collapse cluster could become as large as 4. This is 
consistent with the fact that stellar density profile of 
M15 obtained by HST does not show flat core down 
to 2" (Guhathakurta et al. 1996). Note, however, 
that the core radius obtained from the light profile 
is somewhat larger than the core radius of the mass; 
more importantly, real clusters are tidally limited and 
therefore are not described accurately by these mod- 
els. 
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5. GRAVOTHERMAL OSCILLATIONS 

The inner regions of a slowly evolving postcol- 
lapse cluster are nearly isothermal. Their slow ex- 
pansion on the timescale t r h is almost negligible on 
the local timescale t rc . Thus these regions resem- 
ble equilibrium isothermal spheres, and therefore they 
are subject to the gravothermal instability studied 
by Antonov (1962) and Lynden-Bell & Wood (1968). 
This instability of postcollapse clusters was first found 
by Sugimoto & Bettwieser (1983) and Bettwieser & 
Sugimoto (1984), and later investigated by a number 
of scientists. 

5.1. The Instability Parameter e 

Single-component, isolated clusters in postcollapse 
form a one-parameter family when they are domi- 
nated by three-body binaries. That parameter can 
be taken to be N, the number of stars. In this case, 
gravothermal instability occurs for TV <^ N cr u ~ 7000 
(Goodman 1987). The introduction of multiple com- 
ponents brings additional dimensionless parameters, 
which affect the energy-generation rate in the core 
and perhaps also the relaxation rate near r% . We may 
therefore expect that the instability should depend on 
parameters other than N . 

Goodman (1993) suggested that the quantity 



_ Efot/trh 
Ec /trc 



(24) 



should describe the degree of stability universally (re- 
gardless of the presence of mass spectrum) , where 



2tt 



Pcrlvl 



(25) 



is the energy of the core, and t r h and t rc are half mass 
and core relaxation times, respectively. The motiva- 
tion for this idea is that the core luminosity L c is 
stabilizing: equations ( fl2"| ) and (|l4]) indicate that L c 
will increase as the core shrinks and decrease as it 
expands. The stabilizing influence will be ineffective 
if e <C 1, however, because then the "equilibrium" 
luminosity of the core (oc E to t/t r h) is very small com- 
pared to the rate at which heat can be removed from 
the core if isothermal conditions should break down 
(oc E c /t rc ). 

Our two-component models provide a test of these 
ideas since, as Goodman (1993) argued, e depends 
both on N and on ma/mi. We have performed an- 
other set of runs (Group B) whose parameters were 



chosen to test the analytical predictions given below. 
These simulations include heating by three-body bi- 
naries only. The initial conditions were Plummer 
models with 3 different mass ratios (m 2 /mi = 1.5, 
2, 3) and 4 different total numbers (N = 3 x 10 4 , 1 x 



10 5 ,3 x 10 5 ,1 x 10 6 ). We have fixed N 1 /N 2 



100 



except for four supplementary runs. The parameters 
of the Group B runs are shown in Table ||. 

As in Table | and Figure |, 6 out of 12 runs showed 
gravothermal oscillations in the postcollapse expan- 
sion phase. Figure [?] shows that e itself oscillates. An 
"equilibrium" value for e can be obtained by adopt- 
ing integration time steps larger than the typical os- 
cillation period; with the implicit time integration 
schemes that we and others use, large time steps sup- 
press the gravothermal oscillations. The equilibrium e 
is almost constant during the whole instability period 
(Figure 0) and can be used as a representative value 
for the postcollapse phase. These representative e's 
are plotted in Figure || against a combination of clus- 
ter parameters that will be explained in the following 
subsection. There is a clear boundary of stability at 
e crit 0.01, close to the value 0.013 found by Good- 
man (1993) for single-component clusters. 

It is instructive to compare certain rows in Table 
2. Runs go04, go05, and go06 have the same num- 
ber of stars in the heavier component but increasing 
values of 7712/ 'mi, e, and r c /rh- The first run shows 
oscillations, and the latter two do not. Therefore, the 
suggestion of Murphy, Cohn, & Hut (1990) that sta- 
bility depends on the number of heavy stars is not 
borne out by these runs. Note by the way that the 
number of heavy stars is only 1000, much less than 
the critical value for instability in single-mass clusters 
(7000). It is also interesting to compare the unstable 
run go04 with run go09: the latter has three times as 
many stars in each component but is stable, appar- 
ently because of its larger 7712 / mi . 

5.2. Dependence of e on Cluster Parameters 

In their study on the stability of clusters with 
three-body binaries and a broad mass spectrum, Mur- 
phy et al. (1990) found that stability persists to much 
larger N than in single-mass clusters if the mass func- 
tion is steep. They suggested that stability depends 
mainly on the total number of the heaviest stars. We 
believe that e < 0.01 is a more reliable criterion for 
the onset of instability. 

We now determine e MS H function of cluster param- 
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Table 4 

Parameters and Results of Simulation Group B 



Run 


r»2 
mi 


Ni 
N 2 


N 


e 


r e 
r h 






Oscillation? 


goOl 


1.5 


100 


3 x 10 4 


1.09 x 10~ 2 


4.09 x 


10 


-3 


n 


go02 


2.0 


100 


3 x 10 4 


2.40 x 10~ 2 


7.03 x 


10 


-3 


n 


go03 


3.0 


100 


3 x 10 4 


5.96 x 10~ 2 


1.31 x 


10 


-2 


n 


go04 


1.5 


100 


10 5 


4.97 x 10~ 3 


1.69 x 


10 


-3 


y 


go05 


2.0 


100 


10 5 


1.08 x 10~ 2 


2.94 x 


10 


-3 


11 


go06 


3.0 


100 


10 5 


2.53 x 10~ 2 


5.38 x 


10 


-3 


11 


go07 


1.5 


100 


3 x 10 5 


2.13 x lO" 3 


6.47 x 


10 


-4 


y 


go08 


2.0 


100 


3 x 10 5 


4.61 x 10~ 3 


1.13 x 


10 


-3 


y 


go09 


3.0 


100 


3 x 10 5 


1.12 x 10~ 2 


2.33 x 


10 


-3 


n 


golO 


1.5 


100 


10 6 


7.73 x 10~ 4 


2.09 x 


10 


-4 


y 


goll 


2.0 


100 


10 6 


1.80 x 10~ 3 


4.25 x 


10 


-4 


y 


gol2 


3.0 


100 


10 6 


4.08 x 10~ 3 


8.23 x 


10 


-4 


y 








Supplementary Runs 








go05a 


2.0 


30 


10 5 


0.97 x 10~ 2 


2.70 x 


10 


-3 


y 


go05b 


2.0 


300 


10 5 


0.90 x 10~ 2 


2.57 x 


10 


-3 


n 


go08a 


2.0 


30 


3 x 10 5 


3.98 x 10~ 3 


1.05 x 


10 


-3 


y 


go08b 


2.0 


300 


3 x 10 5 


4.15 x 10~ 3 


1.05 x 


10 


-3 


y 



Note. — e and r c /r^ values are equilibrium values. 
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Fig. 6. — Oscillation-suppressed e values of runs 
in Group B. Filled circles are for runs that showed 
gravothermal oscillation and open circles are for those 
that did not. There is a boundary near e ~ 0.01 below 
which oscillations take place. The diagonal straight 
line represents the simple scaling relation in equation 



d29|). 



eters. We derive simple scaling relations by analytic 
arguments and compare the results with our numeri- 
cal results. 

First we reproduce the prediction by Goodman 
(1993). From equation (]2l|), one obtains 



oc 

rhj \mi 



^ 7 /6 



(26) 



Now, from equations (|ll|) , ([jj]) , and (|l5|) , the ratio of 
energies in the cluster and the core is 



E, 



M z4 r h v% 
oc — — — oc 



E c " M c v 2 c r c vt 



(27) 



where vf n is the velocity dispersion of the whole clus- 
ter. Similarly, the ratio of half-mass to core relaxation 
time is 



rc M, 



1/2 3/2 



x-j — 28 

rr v r m 



where m is the mean mass in the cluster, and m c is 
the mean mass in the core. With rh c /fh w m^jrax 



Fig. 7. — Comparison of oscillation-suppressed 
(dashed lines) and unsuppressed (solid lines) evolu- 
tion of central density (upper lines) and e (lower 
lines) of run gol2. The oscillation-suppressed val- 
ues are obtained by appropriately enlarging integra- 
tion timesteps for the Fokker- Planck equation. The 
oscillation-suppressed values for e are good (and al- 
most constant) representatives of e during oscillation. 
p c is in units of 10 5 Af Q pc~ 3 and t in code units 

(— ^rh.i)- 

and ~ 7712/mi, one finally obtains 



e oc 



oc 



>_c_ 

rh 
m 2 



m 2 
mi 



1/2 



5/3 



— ] N-V\ 
mi 



(29) 



The exponents of mass ratio 3/2 and 2 in equation 
(18) of Goodman (1993) should be corrected to 7/6 
and 5/3 as in equations ( |26| ) and (|29|). 

We have compared the above scaling relations with 
our numerical results from runs in Group B. The equi- 
librium e values of our runs are plotted in Figure || 
against the righthand side of equation (|2^) . The data 
points are well aligned and their slope is 1.20 ± 0.05, 
slightly higher than the theoretical value of unity. 

This discrepancy is traceable to the deviation of 
Tc/fh from the predicted scaling (p6|). In §|[ we re- 
marked that our numerical results show a slightly 
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higher slope than expected for r c (1.07 ± 0.01 ver- 
sus 1) and a lower slope for (0.86 ± 0.06 versus 1). 
Thus the slope for r c /r^ should be higher than pre- 
dicted by our analytic arguments by appproximately 
0.21 for the runs in Group A. The correlation between 
the two sides of equation ( p6| ) for Group B is shown in 
Figure [|. The slope of these data points is 1.26±0.04, 
which is comparable to the one for e. We find that 
one of the causes of this poor r c /r^ approximation is 
the assumption that the proportionality constant in 
equation ( |l6| ) is independent of cluster parameters. 
However, we find that the ratio (rh/t r h)/fh of our 
runs ranges from 10 to 15. Runs with larger mijm\ 
show larger {rh/t r h)/rh values and this mijm\ de- 
pendence is more pronounced for runs with larger N. 
If one lets 

Ar h: (30) 



trh 



where the coefficient A depends on cluster parameters 
such as 777,2/TOi and N, one finds 



r c /r h ocA 8 / 9 . 



(31) 



Then a 23 % variation in A would give a 20 % residual 
in the slope of the correlation between the two sides 
of equation (|l^). This is approximately what one sees 
in Figure (8) , which is based on m 2 /mi ranging from 
1.5 to 3. Although equation (pjl) has been used widely 
without consideration of its dependence on cluster pa- 
rameters, we conclude that the coefficient in this re- 
lation varies with the mass function. 

We have made four more runs in addition to Group 
B to test the dependence of e on Ni/N 2 , which is not 
apparent in equation (p9|): two runs with the same 
parameters as run go05 except Ni/N 2 = 30 and 300 
(go05a and go05b, respectively), and two runs as run 
go08 except N x /N 2 = 30 and 300 (go08a and go08b, 
respectively). Equilibrium e values of these runs are 
0.97 x 10" 2 for run go05a, 0.90 x 10~ 2 for run go05b, 
3.98 x 10~ 3 for run go08a, and 4.15 x 10~ 3 for run 
go08b. These values are all within only 15 % dif- 
ferences from their comparison runs, go05 and go08, 
indicating that e is independent of N1/N2 as expected 
in the above energy balance analysis. However, while 
gravothermal oscillations are observed in run go08a 
and not in runs go05a and go05b as expected, it is 
not observed in run go08b. Clearly the criterion of 
e ;> 0.01 does not appear to be exact if N 2 is too 
small. 

All of our previous analyses assume that there are 
enough stars in the heavy component so that the core 




-4 -3 -2 

log (m/m,) 7 ^- 2 / 3 

Fig. 8. — Oscillation-suppressed r c /rh values of our 
runs in Group B. Filled circles are for runs that 
showed gravothermal oscillation and open circles are 
for those that did not. The straight diagonal line rep- 
resents the theoretical scaling relation. 



is dominated by them. This requires that the ratio of 
core mass to total mass, 



M c /M w 3.3 x 



m 2 
mi 



1/6 



iv- 2 / 3 t n ° 



(15 



(32) 



(see eq. |22e|) should be smaller than M 2 /M. In or- 
der to confine most of the heavy stars to a region 
much smaller than r^, we also require m 2 jm\ > 3/2: 
p 2 (r) oc pi(r) m2 / mi in equipartion; p\ oc r~ 2 in 
the regions well outside the core where the lighter 
component dominates the potential; and we require 
P2{t) to fall more steeply than r~ 3 in those regions 
so that most of the heavies are at smaller radii. 
All of our runs satisfy these requirements, except 
where m 2 /m\ = 1.5 so that the second condition is 
marginally violated. 

6. SUMMARY 

We have investigated the evolution of isolated two- 
component clusters with Plummer-model initial con- 
ditions. We have included the main effects of both 
three-body and tidal-capture binaries by adding a 
heating term to the Fokker-Planck equation. 
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In agreement with previous investigators, we find 
that core collapse is hastened by the presence of heavy 
remnant stars. 

In the postcollapse phase, we find that heating by 
three-body binaries exceeds that by tidal binaries at 
least for clusters with M < 3 x 10 5 M Q , ru,i > 5pc, 
to 2 /toi > 2, and N x /N 2 < 300 when m 2 = 1.4 M . 
When three-body binary heating does dominate, the 
expansion of the postcollapse cluster is self-similar. 
Scaling laws for cluster parameters including the cen- 
tral density, velocity dispersion, core radius, and half- 
mass radius have been derived from simple consider- 
ations of energy balance, and these scalings generally 
agree well with our numerical results, which however 
also provide the numerical coefficients in the scaling 
laws. We related the postcollapse evolution of these 
cluster parameters to N, M and m\jm 2 . 

We have studied the gravothermal oscillation phe- 
nomenon using our two-component models. We have 
confirmed that the parameter e = {E to t/t r h)/(E c /t r h) 
predicts the occurence of gravothermal oscillations in 
the presence of this simplest of nontrivival mass func- 
tions. The scaling law for e with respect to mi/m^ 
and N is derived in the limit of small N 2 /N\ and com- 
pared with our numerical results. Generally speaking, 
clusters with a steeper mass function are less suscep- 
tible to gravothermal instability. This is in qualita- 
tive agreement with an earlier suggestion by Murphy 
et al. (1990), but whereas these authors suggested 
that stability depends on the number of heavy stars, 
we conclude that it is independent of the number of 
heavies but depends jointly on the total number of 
stars and on the ratio of the individual stellar masses 
in the two components. 

These conclusions need to be tested against models 
with more mass components. Our preliminary stud- 
ies indicate that the evolution of multi-mass clusters 
is generally similar to two-component models for an 
appropriate choice of model parameters (mostly suit- 
able m 2 /mi). In order to compare with real clusters, 
we need to take into account a host of complicating 
physical effects, most importantly an external tidal 
field. The results of these studies will be reported 
elsewhere. 
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